Cargamos los paquetes necesarios
library(dplyr)
library(sf)
library(gstat)
library(geoR)
library(spdep)
library(DataExplorer)
library(psych)
library(ggcorrplot)
library(biogeo)
library(ggplot2)
library(leaflet)
Cargamos los datasets
estaciones <- read.csv("data_smn/preprocessed/estaciones_smn_v2.csv")
horarios <- read.csv("data_smn/preprocessed/datohorario20210420.csv")
Observamos como esta compuesto el dataset
glimpse(estaciones)
Rows: 123
Columns: 9
$ NOMBRE <chr> "BASE BELGRANO II", "BASE CARLINI (EX JUBANY)", "BASE ESPERANZA", "BASE MARAMBIO", "BASE O…
$ PROVINCIA <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "BUENOS AIRE…
$ LATITUD_GRADOS <int> 77, 62, 63, 64, 60, 68, 36, 38, 37, 36, 34, 38, 37, 36, 34, 34, 34, 34, 36, 37, 34, 34, 34…
$ LATITUD_MINUTOS <int> 52, 14, 24, 14, 45, 8, 50, 44, 43, 12, 32, 1, 26, 21, 36, 49, 33, 58, 2, 56, 33, 41, 40, 2…
$ LONGITUD_GRADOS <int> 34, 58, 56, 56, 44, 67, 59, 62, 59, 61, 58, 61, 61, 57, 58, 58, 60, 57, 59, 57, 58, 58, 58…
$ LONGITUD_MINUTOS <int> 34, 38, 59, 43, 43, 8, 53, 10, 47, 4, 40, 20, 53, 44, 36, 32, 55, 54, 8, 35, 49, 44, 38, 5…
$ ALTURA <int> 256, 11, 24, 198, 12, 7, 147, 83, 207, 94, 26, 247, 233, 9, 12, 20, 81, 23, 36, 21, 32, 36…
$ NRO <int> 89034, 89053, 88963, 89055, 88968, 89066, 87641, 87750, 87649, 87640, 87570, 87683, 87637,…
$ NroOACI <chr> "SAYB", "SAYJ", "SAYE", "SAWB", "SAYO", "SAYS", "SAZA", "SAZB", "SAZJ", "SAZI", "SADO", "S…
print("#--------------------#")
[1] "#--------------------#"
glimpse(horarios)
Rows: 2,041
Columns: 8
$ FECHA <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, …
$ HORA <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4,…
$ TEMP <dbl> 20.7, 19.8, 19.5, 19.2, 19.2, 18.9, 19.0, 18.8, 19.0, 19.1, 19.9, 20.3, 21.8, 22.9, 23.2, 23.7, 23.2…
$ HUM <int> 70, 79, 79, 81, 81, 84, 84, 84, 84, 84, 81, 82, 68, 72, 68, 66, 71, 68, 75, 72, 73, 68, 68, 68, 90, …
$ PNM <dbl> 1020.8, 1020.9, 1020.9, 1020.4, 1020.0, 1020.3, 1020.7, 1021.1, 1021.6, 1022.2, 1022.4, 1021.9, 1021…
$ DD <int> 70, 70, 70, 70, 70, 50, 70, 50, 50, 50, 50, 50, 50, 70, 70, 70, 70, 70, 90, 90, 90, 90, 70, 70, 50, …
$ FF <int> 28, 28, 26, 26, 22, 22, 17, 17, 15, 19, 20, 24, 20, 19, 20, 17, 15, 15, 19, 20, 22, 22, 31, 35, 17, …
$ NOMBRE <chr> "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AERO…
Observamos que la variable que representa la velocidad del viento en km/h (FF) es de tipo Char, por lo cual debemos convertila en int.
horarios$FF <- as.numeric(horarios$FF)
Observamos el resumen estadistico de las variables de potencial interes
print("altura")
[1] "altura"
summary(estaciones$ALTURA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.0 53.5 145.0 331.5 455.0 3459.0
print("hora")
[1] "hora"
summary(horarios$HORA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 7.00 12.00 12.16 18.00 23.00
print("temperatura")
[1] "temperatura"
summary(horarios$TEMP)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-20.00 16.10 20.20 19.09 23.80 33.60 1
print("humedad")
[1] "humedad"
summary(horarios$HUM)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
18.00 59.00 72.00 71.18 85.00 100.00 4
print("presion atmosferica")
[1] "presion atmosferica"
summary(horarios$PNM)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
960.8 1011.0 1015.1 1075.3 1019.2 3133.0 23
print("direccion del viento")
[1] "direccion del viento"
summary(horarios$DD)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 50.0 70.0 118.9 160.0 990.0 1
print("velocidad del viento")
[1] "velocidad del viento"
summary(horarios$FF)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 9.00 15.00 16.06 22.00 107.00 1
Se puede observar la presencia de observaciones NA’s, por lo cual procedemos a eliminarlas
which(is.na(horarios), arr.ind=TRUE)
row col
[1,] 771 3
[2,] 587 4
[3,] 771 4
[4,] 1876 4
[5,] 1877 4
[6,] 237 5
[7,] 771 5
[8,] 973 5
[9,] 974 5
[10,] 975 5
[11,] 976 5
[12,] 977 5
[13,] 978 5
[14,] 979 5
[15,] 980 5
[16,] 1392 5
[17,] 1393 5
[18,] 1394 5
[19,] 1395 5
[20,] 1396 5
[21,] 1397 5
[22,] 1398 5
[23,] 1399 5
[24,] 1400 5
[25,] 1401 5
[26,] 1402 5
[27,] 1403 5
[28,] 1441 5
[29,] 771 6
[30,] 771 7
horarios[771,]
Como toda la fila para la rioja aero es NA, lla eliminamos
horarios <- horarios[-c(771), ]
rownames(horarios) <- 1:nrow(horarios)
horarios[587,]
Completamos la columna humedad con el promedio
# FORMOSA AERO
formosa <- horarios[horarios$NOMBRE == "FORMOSA AERO",]
which(is.na(formosa), arr.ind=TRUE)
row col
587 4 4
formosa <- na.omit(formosa)
horarios[587,4] = mean(formosa$HUM)
horarios[587,]
remove(formosa)
which(is.na(horarios), arr.ind=TRUE)
row col
1875 1875 4
1876 1876 4
237 237 5
972 972 5
973 973 5
974 974 5
975 975 5
976 976 5
977 977 5
978 978 5
979 979 5
1391 1391 5
1392 1392 5
1393 1393 5
1394 1394 5
1395 1395 5
1396 1396 5
1397 1397 5
1398 1398 5
1399 1399 5
1400 1400 5
1401 1401 5
1402 1402 5
1440 1440 5
horarios[c(1875,1876),]
Completamos con el promedio de humedad calculado para tucuman aero
# TUCUMAN AERO
tucuman_aero <- horarios[horarios$NOMBRE == "TUCUMAN AERO",]
which(is.na(tucuman_aero), arr.ind=TRUE)
row col
1875 1 4
1876 2 4
tucuman_aero <- na.omit(tucuman_aero)
horarios[1875,4] = mean(tucuman_aero$HUM)
horarios[1876,4] = mean(tucuman_aero$HUM)
horarios[c(1875,1876),]
remove(tucuman_aero)
horarios[237,]
# CATAMARCA AERO
catamarca_aero <- horarios[horarios$NOMBRE == "CATAMARCA AERO",]
which(is.na(catamarca_aero), arr.ind=TRUE)
row col
237 2 5
catamarca_aero <- na.omit(catamarca_aero)
horarios[237,5] = mean(catamarca_aero$PNM)
horarios[237,]
remove(catamarca_aero)
horarios[c(972,973,974,975,976,977,978,979),]
# MERCEDES AERO
mercedes_aero <- horarios[horarios$NOMBRE == "MERCEDES AERO",]
which(is.na(mercedes_aero), arr.ind=TRUE)
row col
972 1 5
973 2 5
974 3 5
975 4 5
976 5 5
977 6 5
978 7 5
979 8 5
remove(mercedes_aero)
Vamos a completar los valores faltantes para la presion atmosferica con el promedio de la provincia de corrientes a la que pertenece Mercedes aero
corrientes_aero <- horarios[(horarios$NOMBRE == "CORRIENTES AERO") | (horarios$NOMBRE == "ITUZAINGO") | (horarios$NOMBRE == "MONTE CASEROS AERO") | (horarios$NOMBRE == "PASO DE LOS LIBRES AERO"),]
which(is.na(corrientes_aero), arr.ind=TRUE)
row col
horarios[972,5] = mean(corrientes_aero$PNM)
horarios[973,5] = mean(corrientes_aero$PNM)
horarios[974,5] = mean(corrientes_aero$PNM)
horarios[975,5] = mean(corrientes_aero$PNM)
horarios[976,5] = mean(corrientes_aero$PNM)
horarios[977,5] = mean(corrientes_aero$PNM)
horarios[978,5] = mean(corrientes_aero$PNM)
horarios[979,5] = mean(corrientes_aero$PNM)
horarios[c(972,973,974,975,976,977,978,979),]
remove(corrientes_aero)
horarios[c(1391,1392,1394,1395,1396,1397,1398,1399,1400,1401,1402,1403),]
# RIO GALLEGOS AERO
gallegos_aero <- horarios[horarios$NOMBRE == "RIO GALLEGOS AERO",]
which(is.na(gallegos_aero), arr.ind=TRUE)
row col
1391 1 5
1392 2 5
1393 3 5
1394 4 5
1395 5 5
1396 6 5
1397 7 5
1398 8 5
1399 9 5
1400 10 5
1401 11 5
1402 12 5
gallegos_aero <- na.omit(gallegos_aero)
horarios[1391,5] = mean(gallegos_aero$PNM)
horarios[1392,5] = mean(gallegos_aero$PNM)
horarios[1393,5] = mean(gallegos_aero$PNM)
horarios[1394,5] = mean(gallegos_aero$PNM)
horarios[1395,5] = mean(gallegos_aero$PNM)
horarios[1396,5] = mean(gallegos_aero$PNM)
horarios[1397,5] = mean(gallegos_aero$PNM)
horarios[1398,5] = mean(gallegos_aero$PNM)
horarios[1399,5] = mean(gallegos_aero$PNM)
horarios[1400,5] = mean(gallegos_aero$PNM)
horarios[1401,5] = mean(gallegos_aero$PNM)
horarios[1402,5] = mean(gallegos_aero$PNM)
horarios[1403,5] = mean(gallegos_aero$PNM)
horarios[c(1391,1392,1394,1395,1396,1397,1398,1399,1400,1401,1402,1403),]
remove(gallegos_aero)
which(is.na(horarios), arr.ind=TRUE)
row col
1440 1440 5
horarios[1440,]
# RIO GRANDE
rio_grande_aero <- horarios[horarios$NOMBRE == "RIO GRANDE B.A.",]
which(is.na(rio_grande_aero), arr.ind=TRUE)
row col
1440 26 5
rio_grande_aero <- na.omit(rio_grande_aero)
horarios[1440,5] = mean(rio_grande_aero$PNM)
horarios[1440,]
remove(rio_grande_aero)
Validamos que efectivamente se hayan arreglado
print("temperatura")
[1] "temperatura"
summary(horarios$TEMP)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-20.00 16.10 20.20 19.09 23.80 33.60
print("humedad")
[1] "humedad"
summary(horarios$HUM)
Min. 1st Qu. Median Mean 3rd Qu. Max.
18.00 59.00 72.00 71.19 85.00 100.00
print("presion atmosferica")
[1] "presion atmosferica"
summary(horarios$PNM)
Min. 1st Qu. Median Mean 3rd Qu. Max.
960.8 1011.1 1015.1 1074.7 1019.1 3133.0
print("direccion del viento")
[1] "direccion del viento"
summary(horarios$DD)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 50.0 70.0 118.9 160.0 990.0
print("velocidad del viento")
[1] "velocidad del viento"
summary(horarios$FF)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 9.00 15.00 16.06 22.00 107.00
Perfecto, ya no tenemos valores nulos en nuestro dataset.
A continuación, visualizaremos las distribuciones de las variables numericas de interés
hist(estaciones$ALTURA, main = "Histograma de la altura", xlab = "Altura")
hist(horarios$HORA, main = "Histograma de horarios", xlab = "Horarios")
hist(horarios$TEMP, main = "Histograma de temperatura", xlab = "Temperatura")
hist(horarios$HUM, main = "Histograma de humedad", xlab = "Humedad")
hist(horarios$PNM, main = "Histograma de presion atmosferica", xlab = "Presion Atmosferica")
hist(horarios$DD, main = "Histograma de la direccion del viento", xlab = "Direccion del viento")
hist(horarios$FF, main = "Histograma de la velocidad del viento", xlab = "Velocidad del viento")
Creamos unos boxplot para visualizar la distribucionde datos que potencialmente nos interecen para proceder con el analisis
boxplot(estaciones$ALTURA, main = "Boxplot de la altura", xlab = "Altura")
boxplot(horarios$HORA, main = "Boxplot de horarios", xlab = "Horarios")
boxplot(horarios$TEMP, main = "Boxplot de temperatura", xlab = "Temperatura")
boxplot(horarios$HUM, main = "Boxplot de humedad", xlab = "Humedad")
boxplot(horarios$PNM, main = "Boxplot de presion atmosferica", xlab = "Presion Atmosferica")
boxplot(horarios$DD, main = "Boxplot de la direccion del viento", xlab = "Direccion del viento")
boxplot(horarios$FF, main = "Boxplot de la velocidad del viento", xlab = "Velocidad del viento")
Los boxplots anteriores ponen en evidencia la existencia de outliers. ¿Pero son estos realmente outliers, o pertenecen a observaciones en lugares muy remotos? Esto lo analizaremos luego, al momento de graficar las estaciones en el mapa de Argentina.
Ahora, veamos que tan correlacionadas estan estas variables.
corr <- cor(horarios[, c(2,3,4,5,6,7)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
lab=TRUE,
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726"))
Las variables que mas correlacionan con la velocidad del viento son HUMEDAD (negativamente) y HORA (positivamente). Tambien vemos que HORA y TEMPERATURA correlacionan negativamentecon HUMEDAD. Por ultimo, se observa que HORA y TEMPERATURA correlacionan positivamente
Analizamos ahora la simetria de la variable que representa la velocidad del viento ya que es la que mas nos interesa en este estudio.
skew(horarios$FF)
[1] 1.588087
kurtosi(horarios$FF)
[1] 7.188948
La medida de asimetria y kurtosi terminan de validar lo que observamos en el histograma. La variable FF es asimetrica a derecha y tiene una mayor concentracion de valores muy cerca de la media de la distribución y muy lejos de la cola de la distribucion.
Un detalle no menor del dataset de estaciones es que las latitudes y longitudes estan expresadads en grados y minutos. Para poder trabajar con ellas, necesitamos que esten expresadas en valores decimales. Por eso, en el siguiente bloque de código vamos a usar la funcion dms2dd para hacer esta conversión.
# creamos dos vectores vacios
latitud <- c()
longitud <- c()
# iteramos por cada fila del dataset de estaciones y hacemos la convesion de latitud y longitud
for(i in 1:nrow(estaciones)) {
latitud[i] <- dms2dd(dd = estaciones[i, "LATITUD_GRADOS"], mm = estaciones[i, "LATITUD_MINUTOS"], ss = 0, ns = "S")
longitud[i] <- dms2dd(dd = estaciones[i, "LONGITUD_GRADOS"], mm = estaciones[i, "LONGITUD_MINUTOS"], ss = 0, ns = "S")
}
# asignamos a latitud y longitud los valores convertidos
estaciones['LATITUD'] <- latitud
estaciones['LONGITUD'] <- longitud
Antes de unir estaciones, se removeran las columnas que no sean relevantes para este analisis. Las mismas son NRO y NroOACI, LATITUD_GRADOS, LATITUD_MINUTOS, LONGITUD_GRADOS, LONGITUD_MINUTOS. Asi como tambien se removera la variable fecha, ya que estos datos pertenecen al 20/04/2021
estaciones <- estaciones[c(1,2,7,10,11)]
horarios <- horarios[,c(3,4,5,6,7,8)]
summary(estaciones)
NOMBRE PROVINCIA ALTURA LATITUD LONGITUD
Length:123 Length:123 Min. : 3.0 Min. :-77.87 Min. :-72.05
Class :character Class :character 1st Qu.: 53.5 1st Qu.:-38.17 1st Qu.:-66.47
Mode :character Mode :character Median : 145.0 Median :-34.53 Median :-63.37
Mean : 331.5 Mean :-35.96 Mean :-62.82
3rd Qu.: 455.0 3rd Qu.:-30.25 3rd Qu.:-58.93
Max. :3459.0 Max. :-22.10 Max. :-34.57
summary(horarios$TEMP)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-20.00 16.10 20.20 19.09 23.80 33.60
Se procede a unificar las dos tablas usando la variable NOMBRE como punto para combinar los datasets
data <- inner_join(estaciones, horarios, by = c("NOMBRE" = "NOMBRE"))
glimpse(data)
Rows: 2,040
Columns: 10
$ NOMBRE <chr> "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II…
$ PROVINCIA <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTAR…
$ ALTURA <int> 256, 256, 256, 256, 256, 256, 256, 256, 11, 11, 11, 11, 11, 11, 11, 11, 24, 24, 24, 24, 24, 24, 2…
$ LATITUD <dbl> -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -62.23333…
$ LONGITUD <dbl> -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -58.63333…
$ TEMP <dbl> -20.0, -16.2, -16.9, -13.0, -15.1, -16.6, -15.3, -16.8, -0.1, -0.5, -0.9, -0.5, 0.3, 1.1, 1.3, 2.…
$ HUM <dbl> 70, 72, 80, 73, 86, 66, 77, 85, 88, 82, 77, 75, 84, 92, 99, 91, 43, 49, 67, 81, 67, 62, 77, 83, 7…
$ PNM <dbl> 975.7, 970.4, 966.7, 964.2, 961.4, 963.0, 964.8, 967.1, 992.9, 993.0, 992.9, 994.6, 995.3, 995.6,…
$ DD <int> 140, 140, 110, 140, 160, 160, 90, 110, 230, 230, 270, 290, 290, 320, 320, 320, 160, 200, 50, 200,…
$ FF <dbl> 13, 30, 50, 52, 74, 61, 44, 13, 15, 15, 26, 28, 41, 37, 50, 74, 20, 7, 13, 7, 11, 46, 56, 74, 19,…
summary(data$NOMBRE)
Length Class Mode
2040 character character
summary(data$PROVINCIA)
Length Class Mode
2040 character character
summary(data$ALTURA)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.0 43.0 125.0 309.2 450.0 3459.0
summary(data$LATITUD)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-77.87 -37.63 -34.13 -35.47 -30.27 -22.10
summary(data$LONGITUD)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-72.05 -66.28 -63.02 -62.85 -58.73 -34.57
summary(data$TEMP)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-20.00 16.10 20.20 19.09 23.80 33.60
summary(data$HUM)
Min. 1st Qu. Median Mean 3rd Qu. Max.
18.00 59.00 72.00 71.19 85.00 100.00
summary(data$PNM)
Min. 1st Qu. Median Mean 3rd Qu. Max.
960.8 1011.1 1015.1 1074.7 1019.1 3133.0
summary(data$DD)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 50.0 70.0 118.9 160.0 990.0
summary(data$FF)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 9.00 15.00 16.06 22.00 107.00
Volvemos realizar el grafico de correlacion incluyendo la variable altura.
corr <- cor(data[, c(3,4,5,6,7,8,9,10)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
lab=TRUE,
ggtheme = ggplot2::theme_gray,
colors = c("#6D9EC1", "white", "#E46726"))
Se observa que hay estaciones que tiene las observaciones en 0 para la variable FF y DD. Consideramos esto como un error en el instrumento de medicion, por lo cual vamos a eliminar a esa estacion del analisis.
data = data[(data$FF != 0) & (data$DD != 0),]
rownames(data) <- 1:nrow(data)
Agrupamos los datos por nombre calculando el promedio y desvio del viento
data_agg = data %>%
group_by(NOMBRE) %>%
summarise(MEAN_VIENTO_KMH = mean(FF),
SD_VIENTO_KMH = sd(FF),
LONGITUD = unique(LONGITUD),
LATITUD = unique(LATITUD),
.groups = "keep")
Vemos que en el dataset resultante nos quedan 98 observaciones que coinciden con la cantidad de estaciones meteorologicas originales
Convertimos data_agg a data.frame ya que necesitamos este tipo de dato para poder trabajar
data_agg = data.frame(data_agg)
Transformamos df_data_agg en un archivo geográfico utilizando el código de proyección mercator
data_agg_sf = st_as_sf(data_agg, coords = c("LONGITUD", "LATITUD"), crs = 4326)
Validamos la clase del nuevo dataframe
class(data_agg_sf)
[1] "sf" "data.frame"
En el siguiente gráfico de la republica argentina se observan en color azul las estaciones meteorológicas donde se realizaron las mediciones de la variables que estan presentes en el dataset
Queremos que en el mapa se vea como etiqueta el nombre de la base meteorologica. Para eso aplicamos la siguiente funcion
labs <- lapply(seq(nrow(data_agg_sf)), function(i) {
paste0( '<p>', data_agg_sf[i, "NOMBRE"], '<p>', '<p>',data_agg_sf[i, "MEAN_VIENTO_KMH"],'</p>' ) })
Realizamos el grafico interactivo de las estaciones meteorologicas graficadas sobre un mapa de Argentina.
leaflet() %>%
addTiles() %>%
addCircles(data = data_agg_sf, weight = 3, label = lapply(labs, htmltools::HTML))
NA
En el mapa observamos que hay puntos muy distantes de la Argentina continental. Dado el proposito de este estudio, el cual es determinar la ubicacion geografica óptima en base a la variable velocidad del viento, decidimos remover estas observaciones ya que no aporta informacion util y ademas agregan ruido a nuestro analisis.
Primero, vamos a borrar las estaciones que no estan en la plataforma continental argentina - Base Carlini - Base San Martin - Base Marambio - Base Esperanza - Base Orcadas
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "BASE CARLINI (EX JUBANY)",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "BASE SAN MARTIN",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "BASE MARAMBIO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "BASE ESPERANZA",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "BASE ORCADAS",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "BASE BELGRANO II",]
Repetimos el plot para validar
labs <- lapply(seq(nrow(data_agg_sf)), function(i) {
paste0( '<p>', data_agg_sf[i, "NOMBRE"], '</p>' ) })
leaflet() %>%
addTiles() %>%
addCircles(data = data_agg_sf, weight = 3, label = lapply(labs, htmltools::HTML))
NA
Observemos el resumen estadistico de las nuevas variables MEAN_VIENTO_KMH y SD_VIENTO_KMH
describe(data_agg_sf$MEAN_VIENTO_KMH)
hist(data_agg_sf$MEAN_VIENTO_KMH)
boxplot(data_agg_sf$MEAN_VIENTO_KMH)
describe(data_agg_sf$SD_VIENTO_KMH)
hist(data_agg_sf$SD_VIENTO_KMH)
boxplot(data_agg_sf$SD_VIENTO_KMH)
Veamos si esta nueva variable MEAN_VIENTO_KMH es normal.
hist(data_agg_sf$MEAN_VIENTO_KMH)
boxplot(data_agg_sf$MEAN_VIENTO_KMH)
qqnorm(data_agg_sf$MEAN_VIENTO_KMH)
qqline(data_agg_sf$MEAN_VIENTO_KMH, col=2)
shapiro.test(data_agg_sf$MEAN_VIENTO_KMH)
Shapiro-Wilk normality test
data: data_agg_sf$MEAN_VIENTO_KMH
W = 0.94243, p-value = 9.648e-05
Claramente la variable MEAN_VIENTO_KMH no es normal. El qqplot pone en evidencia la existencia de colas pesadas. Ademas, al realizar el test de shapiro wilk el p-value obtenido es menor a 0.05, lo cual indica que los datos que tenemos no son normales #TODO: agregar algo mas a esta conclusion?
Ahora, procedemos a analizar la existencia de inliers, y en el caso de encontrarlos, eliminarlos.Para eso, usamos el test de moran. Basicamente lo que estamos testeando es que el promedio del viento este dristribuido de manera aleatoria siguiendo un proceso aleatorio.
knea <- knearneigh(data_agg_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)
moran_test <- moran.test(data_agg_sf$MEAN_VIENTO_KMH, listw)
moran_test
Moran I test under randomisation
data: data_agg_sf$MEAN_VIENTO_KMH
weights: listw
Moran I statistic standard deviate = 3.8355, p-value = 6.265e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.423564794 -0.008849558 0.012710177
geary_test <- geary.test(data_agg_sf$MEAN_VIENTO_KMH, listw)
geary_test
Geary C test under randomisation
data: data_agg_sf$MEAN_VIENTO_KMH
weights: listw
Geary C statistic standard deviate = 4.083, p-value = 2.223e-05
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.46134647 1.00000000 0.01740455
shaphiro_test <- shapiro.test(data_agg_sf$MEAN_VIENTO_KMH)
shaphiro_test
Shapiro-Wilk normality test
data: data_agg_sf$MEAN_VIENTO_KMH
W = 0.94243, p-value = 9.648e-05
moran <- moran.plot(data_agg_sf$MEAN_VIENTO_KMH, listw = listw)
data_agg_sf[10,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -65.76667 ymin: -28.6 xmax: -65.76667 ymax: -28.6
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
16 CATAMARCA AERO 47.5 11.19949 POINT (-65.76667 -28.6)
data_agg_sf[31,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -63.75 ymin: -35.7 xmax: -63.75 ymax: -35.7
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
37 GENERAL PICO AERO 30.25 4.234777 POINT (-63.75 -35.7)
data_agg_sf[79,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -64.23333 ymin: -33.11667 xmax: -64.23333 ymax: -33.11667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
85 RIO CUARTO AERO 32.5 7.607014 POINT (-64.23333 -33.11667)
data_agg_sf[108,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -61.96667 ymin: -33.66667 xmax: -61.96667 ymax: -33.66667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
114 VENADO TUERTO AERO 32 11.09955 POINT (-61.96667 -33.66667)
data_agg_sf[68,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -62.38333 ymin: -37.6 xmax: -62.38333 ymax: -37.6
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
74 PIGUE AERO 38.23077 6.905962 POINT (-62.38333 -37.6)
data_agg_sf[21,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -61.88333 ymin: -37.43333 xmax: -61.88333 ymax: -37.43333
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
27 CORONEL SUAREZ AERO 27.33333 5.278889 POINT (-61.88333 -37.43333)
data_agg_sf[89,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -66.35 ymin: -33.26667 xmax: -66.35 ymax: -33.26667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
95 SAN LUIS AERO 17 7.901568 POINT (-66.35 -33.26667)
data_agg_sf[23,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -57.73333 ymin: -36.35 xmax: -57.73333 ymax: -36.35
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
29 DOLORES AERO 15.83333 6.105355 POINT (-57.73333 -36.35)
data_agg_sf[109,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -65.43333 ymin: -36.21667 xmax: -65.43333 ymax: -36.21667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
115 VICTORICA 14 8.660254 POINT (-65.43333 -36.21667)
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "CATAMARCA AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "GENERAL PICO AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "RIO CUARTO AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "VENADO TUERTO AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "PIGUE AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "CORONEL SUAREZ AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "SAN LUIS AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "DOLORES AERO",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "VICTORICA",]
rownames(data_agg_sf) <- 1:nrow(data_agg_sf)
# Creamos una lista de vecinos
knea <- knearneigh(data_agg_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)
# Hacemos el test de moran
moran_test_v2 <- moran.test(data_agg_sf$MEAN_VIENTO_KMH, listw)
moran_test_v2
Moran I test under randomisation
data: data_agg_sf$MEAN_VIENTO_KMH
weights: listw
Moran I statistic standard deviate = 4.1509, p-value = 1.656e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.483041908 -0.009615385 0.014086783
geary_test_v2 <- geary.test(data_agg_sf$MEAN_VIENTO_KMH, listw)
geary_test_v2
Geary C test under randomisation
data: data_agg_sf$MEAN_VIENTO_KMH
weights: listw
Geary C statistic standard deviate = 4.6029, p-value = 2.083e-06
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.41937891 1.00000000 0.01591176
shaphiro_test_v2 <- shapiro.test(data_agg_sf$MEAN_VIENTO_KMH)
shaphiro_test_v2
Shapiro-Wilk normality test
data: data_agg_sf$MEAN_VIENTO_KMH
W = 0.96597, p-value = 0.008459
moran.plot(data_agg_sf$MEAN_VIENTO_KMH, listw = listw)
qqnorm(data_agg_sf$MEAN_VIENTO_KMH)
qqline(data_agg_sf$MEAN_VIENTO_KMH, col=2)
data_agg_sf[56,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -55.13333 ymin: -27.48333 xmax: -55.13333 ymax: -27.48333
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
56 OBERA 4 0 POINT (-55.13333 -27.48333)
data_agg_sf[96,]
Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -62.73333 ymin: -35.96667 xmax: -62.73333 ymax: -35.96667
Geodetic CRS: WGS 84
NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH geometry
96 TRENQUE LAUQUEN 10.33333 3.05505 POINT (-62.73333 -35.96667)
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "OBERA",]
data_agg_sf = data_agg_sf[data_agg_sf$NOMBRE != "TRENQUE LAUQUEN",]
rownames(data_agg_sf) <- 1:nrow(data_agg_sf)
# Creamos una lista de vecinos
knea <- knearneigh(data_agg_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)
# Hacemos el test de moran
moran_test_v2 <- moran.test(data_agg_sf$MEAN_VIENTO_KMH, listw)
moran_test_v2
Moran I test under randomisation
data: data_agg_sf$MEAN_VIENTO_KMH
weights: listw
Moran I statistic standard deviate = 4.8525, p-value = 6.096e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.573582602 -0.009803922 0.014453793
geary_test_v2 <- geary.test(data_agg_sf$MEAN_VIENTO_KMH, listw)
geary_test_v2
Geary C test under randomisation
data: data_agg_sf$MEAN_VIENTO_KMH
weights: listw
Geary C statistic standard deviate = 5.2198, p-value = 8.958e-08
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.33433216 1.00000000 0.01626356
shaphiro_test_v2 <- shapiro.test(data_agg_sf$MEAN_VIENTO_KMH)
shaphiro_test_v2
Shapiro-Wilk normality test
data: data_agg_sf$MEAN_VIENTO_KMH
W = 0.96491, p-value = 0.007798
moran.plot(data_agg_sf$MEAN_VIENTO_KMH, listw = listw)
qqnorm(data_agg_sf$MEAN_VIENTO_KMH)
qqline(data_agg_sf$MEAN_VIENTO_KMH, col=2)
# CALCULAMOS EL COEFICIENTE DE MORAN Y EL DE GEARY
moran(data_agg_sf$MEAN_VIENTO_KMH, listw, length(listw$weights),Szero(listw),zero.policy = FALSE)
$I
[1] 0.5735826
$K
[1] 2.304703
#TODO: Agregar conclusion aca
plot(data_agg_sf$MEAN_VIENTO_KMH)
v <- variogram(MEAN_VIENTO_KMH~1, data_agg_sf)
plot(v)
vt_exp = fit.variogram(v, vgm(125, "Exp", 30, 5))
vt_exp
plot(v , vt_exp)
vt_mat = fit.variogram(v, vgm(125, "Mat", 30, 5))
plot(v , vt_mat)
vt_exc = fit.variogram(v, vgm(125, "Exc", 30, 5))
No convergence after 200 iterations: try different initial values?
plot(v , vt_exc)
vt_bes = fit.variogram(v, vgm(125, "Bes", 30, 5))
plot(v , vt_bes)
attr(vt_exp, 'SSErr')
[1] 0.2612195
attr(vt_mat, 'SSErr')
[1] 0.2612195
attr(vt_exc, 'SSErr')
[1] 2.621681
attr(vt_bes, 'SSErr')
[1] 0.3501463
Kriging
departamentos <- st_read("data_departamentos/Codgeo_Pais_x_dpto_con_datos/pxdptodatosok.shp")
Reading layer `pxdptodatosok' from data source `/Users/antonellaschiavoni/Antonella/Maestria/EstadisticaEspacial/TPFinal/Datos/trabajo_final/smn/estadistica_espacial/data_departamentos/Codgeo_Pais_x_dpto_con_datos/pxdptodatosok.shp' using driver `ESRI Shapefile'
Simple feature collection with 527 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.02985 ymin: -90 xmax: -25.02314 ymax: -21.74506
Geodetic CRS: WGS 84
departamentos <-departamentos[departamentos$departamen != "Antártida Argentina",]
departamentos <-departamentos[departamentos$departamen != "Islas del Atlántico Sur",]
departamentos <- as_Spatial(departamentos)
grilla <- as.data.frame(spsample(departamentos, type="regular", n=5000))
CRS object has comment, which is lost in output
names(grilla) <- c("X", "Y")
coordinates(grilla) <- c("X", "Y")
plot(grilla)
gridded(grilla) <- TRUE
grid has empty column/rows in dimension 2
fullgrid(grilla) <- TRUE
plot(grilla)
proj4string(grilla) <- proj4string(departamentos)
CRS object has comment, which is lost in output
data_agg_sf <- as_Spatial(data_agg_sf)
proj4string(data_agg_sf) <- proj4string(departamentos)
CRS object has comment, which is lost in outputCRS object has comment, which is lost in output
ko1 <- krige(MEAN_VIENTO_KMH~1, data_agg_sf, grilla, model = vt_exp, nmax=20)
[using ordinary kriging]
spplot(ko1["var1.pred"])
spplot(ko1["var1.var"])
ko2 <- krige(MEAN_VIENTO_KMH~1, data_agg_sf, grilla, model = vt_exp, nmax=50)
[using ordinary kriging]
spplot(ko2["var1.pred"])
spplot(ko2["var1.var"])
ko3 <- krige(MEAN_VIENTO_KMH~1, data_agg_sf, grilla, model = vt_bes, nmax=50)
[using ordinary kriging]
spplot(ko3["var1.pred"])
spplot(ko3["var1.var"])
library(raster)
Attaching package: ‘raster’
The following object is masked from ‘package:dplyr’:
select
r <- raster(ko1)
r.m <- mask(r, departamentos)
library(tmap)
tm_shape(r.m) +
tm_raster(n=10,
palette="Blues",
auto.palette.mapping=FALSE,
title="") +
tm_legend(legend.outside=TRUE)
The argument auto.palette.mapping is deprecated. Please use midpoint for numeric data and stretch.palette for categorical data to control the palette mapping.
r <- raster(ko1, layer="var1.var")
r.m <- mask(r, departamentos)
tm_shape(r.m) +
tm_raster(n=7,
palette ="Reds",
title="Variance map ") +
tm_legend(legend.outside=TRUE)
r <- sqrt(raster(ko1, layer="var1.var")) * 1.96
r.m <- mask(r, departamentos)
tm_shape(r.m) +
tm_raster(n=7,
palette ="Reds",
title="95% CI map \n(en km/h)") +
tm_legend(legend.outside=TRUE)